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Abstract 

The scientific objectives of the Lisa Technology Package (LTP) experiment on board of the LISA 
Pathfinder mission demand for accurate cahbration and vaUdation of the data analysis tools in 
advance of the mission launch. The level of confidence required in the mission outcomes can be 
reached only by intensive testing the tools on synthetically generated data. A flexible procedure 
allowing the generation of cross-correlated stationary noise time series was set-up. Multichannel 
time series with the desired cross-correlation behavior can be generated once a model for a mul- 
tichannel cross-spectral matrix is provided. The core of the procedure comprises a noise coloring, 
multichannel filter designed via a frequency-by-frequency eigen-decomposition of the model cross- 
spectral matrix and a subsequent fit in the Z-domain. The common problem of initial transients 
in filtered time series is solved with a proper initialization of the filter recursion equations. The 
noise generator performance was tested in a two-dimensional case study of the closed-loop LTP 
dynamics along the two principal degrees-of-freedom. 
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Introduction 

The LTP (LISA Technology Package) experiment is the main scientific payload on the 
European Space Agency mission, LISA Pathfinder. Its goal is to determine and analyse 
all possible sources of disturbance which perturb the free-falling test masses from their 
geodesic motion. The system is composed of two test masses whose position is sensed by 
an interferometer. The spacecraft cannot simultaneously follow both masses, and so the 
trajectory of only one test mass serves as the drag-free reference along the x (measurement) 
axis. To prevent the trajectories of the two masses from diverging in response to any 
differential force, the second test mass is electrostatically actuated to follow the spacecraft. 
In the main science operating mode, the position of the SC relative to the first test mass 
is controlled using micro-Newton thrusters attached to the SC. The position of the second 
test mass is controlled using capacitive actuators surrounding the test mass. The first 
interferometer channel measures the position of the first test mass relative to the spacecraft. 
The second interferometer channel (differential channel) measures the relative displacement 
between the two test masses. 

A set of different experiments, such as measurements of parasitic voltages, test mass 
charging, thermal and magnetic disturbances, completely covers the scheduled 90 days of 
LTP operations; the overall aim of the experiments is to reach the best free-fall quality in a 
step by step procedure in which the result of the previous experiment is used to define the 
detailed configuration of the following experiment. This cascade-like process aims to demon- 
strate the ability to put a test mass into free-fall at a level where any residual acceleration 
is below 3 x 10^^^ m s^^/v^Hz at frequencies around 1 mHz [l-5j. 

Such a demanding accuracy requires a careful calibration of the spectral estimation algo- 
rithms so as to avoid any systematic bias in the estimation of the spectrum of the residual 
acceleration. Due to the limited time duration of the mission, the amount of data available 
will be not enough for a meaningful and robust calibration of the dedicated data analysis 
tools. The natural way to solve the problem is to calibrate and test the tools in advance of 
the mission, by an in-depth analysis of synthetic noise data. The experiment has a total of 
18 measuring channels sensing the movement of the test masses, many of which are coupled 
so that information is contained not only in the individual power-spectra, but also in the 
cross spectral densities between different channels. In order to set-up a reliable test bench 
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for such a system, a robust and flexible multichannel noise generator is required. 

The problem of generating a sequence of random variables having some deflnite statistical 
properties is well examined in literature. Stein and Storer [S] proposed a procedure for which 
the computation of sample values requires the eigen-decomposition of the covariance ma- 
trix of the process. This allows the identiflcation of a transformation matrix that multiplies 
a vector of independent samples to provide a noise series with the desired correlation prop- 
erties. Levin [7] suggested instead to pass white noise through a digital noise coloring fllter 
with a rational transfer function. The problem connected with the initial transient is solved 
with the calculation of K consecutive {K is the order of the noise shaping fllter) output 
values having the same statistical properties as if they were produced by steady-state oper- 
ation of the fllter. An alternative method for fllter initialization was indicated by Kay [8] 
who realized that one has just to specify the initial state for the fllter. The method for the 
calculation of the initial state is based on the Levinson-Durbin algorithm. As an alterna- 
tive, Franklin [U] described a procedure for the simulation of stationary and non- stationary 
Gaussian random processes. The procedure for a non-stationary process is similar to that 
reported in [6| and is based on the Grout factorization of the covariance matrix of the pro- 
cess. The output for the stationary case is, instead, the result of a simulation of a continuous 
system by means of the state space formalism. The initial state is calculated with a linear 
transformation from uncorrelated random noise samples. All the methods in literature deal 
with the generation of a single channel of data with a given correlation function or, analo- 
gously, with a given spectrum. The algorithms proposed by Levin, Kay and Franklin can be 
crudely summarized in three steps; 1) identiflcation of the desired system, 2) initialization 
of the data sequence, 3) generation of the colored noise data sequence from a sequence of 
zero mean, delta correlated random numbers. 

The method proposed in the present paper follows this classical scheme with the relevant 
difference that it is designed to work with multichannel systems (i.e. multiple inputs, mul- 
tiple outputs). In the following, the mathematical basis of the method is developed, and a 
case study is discussed in order to quantitatively assess the reliability of the procedure. 



3 



I. PRINCIPLES OF THE NOISE GENERATION PROCEDURE 



It can be assumed, in complete generality, that the noise to be generated x (t) has a power 
spectral density (PSD) that can be written as: 

S,,{u) = \Hiu)fSo. (1) 

The process x (t) can be thought as the output of a rational continuous filter with transfer 
function H{u), at the input of which is a white, zero-mean noise e{t) with PSD equal to 
^o- The filter H (u) can be written as: 

TV 

H{u) = y^^^, (2) 
with Vh the residue of H (co) in ph [23j. 

Equation (|2| shows that the process x (t) can in turn be considered as: 



N 

X ['] 

h=l 

where: 



it) = J2y'^{t), (3) 



yhit)=r, j eP-'e{t-t')dt'. (4) 



Thus, generating the process x (t) is equivalent to generating the correlated processes 
i/h (t). A discrete process with time-step T can be realized from the recursive equations: 



yt (t + T) = Vh it) t^^^ + e/, (t + T) 



fc=l,...,]V. (5) 

Since the procedure must not diverge, only poles ph with a negative real part (stable 
poles) are considered. The processes eh (t + T) are not independent but it can be verified 
that their cross-correlations is vanishing for time intervals larger than T. Then, indicating 
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with k and m integers values, the cross-correlation between input processes can be written 
as: 



(e. (kT) 6, (mT)) = 6,^^^ [e(^'+^^)^ - l] , (6) 

Pi +Pj 

where the notation () indicates the expectation value operator. If the process e (t) is zero- 
mean Gaussian, then so is ej (kT), and equation ^ is sufficient to determine the statistics. 

The generated series x (kT) exactly represents the sampling of the continuous process 
X (t) with a time-step T, and therefore it also reproduces the aliasing if the filter H {u) has 
a response different from zero at frequencies larger than 1/2T. This is a consequence of the 
well known relation between the spectra of discrete and continuous processes [10] : 



Sxx = (uj + ^ j 



, , 27ik 

1^1 < -jr- (7) 

Here T is the sampling time, S^^ (u) and S^x (i^) are the spectra of the discrete and 
continuous processes respectively. 

From the point of view of the numerical implementation, it is more convenient to start 
from the assumption of a discrete filter. In the case of finite length discrete time series, 
equations ([T]) and ^ can be rewritten as: 



Sxxin) = \Hin)fSo. 

N 

H{Q) = y- (8) 

where Q is the normalized angular frequency ^2 = 27r f/fg, and fg is the sampling fre- 
quency. 

Each element of the partial-fraction expansion in equation ([s]) can be considered as a 
simple autoregressive moving average filter for which the usual recursive relation holds [H] : 



X (n) + aix (n — 1) -|- • • • -|- qnx {n — N) = h^i {n) + ■ ■ ■ + bui iji — M) 
ai = —p and 02, . . . , ajv = 

60 = r and 61, . . . , 6m = 0. (9) 



Here, are the coefficients of the denominator polynomial, bj are the coefficients of the 
numerator polynomial, r and p are residues and poles as written in equation ([s]) and i{n — k) 
is the step k delayed input to the system. Thus the complete noise generation process is 
obtained by: 



N 



x{n) = y^Xfc {n) 



k=l 



Xk in) = pkXk (n-l) + The (n) . (10) 



Each Xk (n) can be generated according to the recursive equation ( 10 ) starting from the 
same white noise series e (n) . 

Such a procedure provides a noise series whose PSD is an accurate replica of the contin- 
uous noise spectrum up to the Nyquist frequency. Aliasing is not reproduced in the discrete 
case. In the rest of the paper, the detailed calculations for the implementation of a discrete 
multichannel procedure are presented. 



II. MULTICHANNEL NOISE GENERATION 



A. Noise Coloring Filter Identification 



A multichannel sequence can be described by the M-dimensional vector: 



(11) 



\yM{t) J 

where yi (t) is the data sequence at the zth channel. If the process is stationary, the 
cross-correlation matrix at a given delay r is defined as [T2] : 



R(t) 



oo 

j yWy^ it + r)dt. 



(12) 



The elements of the matrix R (r) are the cross-correlations between the different elements 
of the multichannel sequence. The symbol f indicates a matrix conjugate transpose. The 
cross-spectral density matrix for the given multichannel process is defined as the Fourier 
transform of the cross-correlation matrix: 
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SM 



R (r) exp (— zwr) dr 



Sii (w) 



Sim (w) 



(13) 



y Smi (t^) ■ • • 'S'j\/A/ (w) J 

A noise coloring multichannel filter is a linear operation which transforms a delta corre- 
lated unitary variance multichannel random sequence (multichannel white noise e (t)) in to 
a noise sequence with the given cross-spectral density matrix. 



N ^ 

Vi (^) = XI / (^) (t - ^) 



(14) 



where hij (r) is the impulse response of the filter between the j'th input and the ith 
output. Assuming that the number of input channels N is the same as the number of 
output channels M, the multichannel coloring filter can be represented by a square matrix. 
The cross-spectral matrix of the output process can be obtained by the combination of the 
cross-spectral matrix of the input and the frequency response of the filter: 



S (w) = H (w) ■ I ■ 



CO 



(15) 



Here I is the unit matrix corresponding to the cross-spectral matrix of the input multi- 
channel white noise process e (t) and H (u) is the frequency response matrix of the multi- 
channel filter. The problem of the generation of a multichannel noise series with the given 
cross-spectral matrix starts from the identification of the noise coloring filter H {u). 

The eigendecomposition of the cross-spectral matrix S (u) is defined as: 



S (cj) = V (oj) ■ S (w) • V 



-1 



(16) 



where V (u) and S [u) are the eigenvector and eigenvalue matrices of the cross-spectral 
matrix S (u). Since S (u) is Hermitian, its eigenvector matrix is unitary, i.e. V (u) V"'' (w) = 



I. Therefore, combining equations (15) and (16), the noise coloring filter can be obtained: 



H (w) = V (w) ■ (tu). (17) 

As S (cj) is a diagonal matrix, a/S (cj) is a diagonal matrix with elements given by the 
square root of the elements of S {u). 



B. System Discretization 

Once the frequency response of the coloring filter H (u) is known, a discrete multichannel 
filter is required for the generation of discrete synthetic noise data series. Discrete filters 
can be estimated by a least square fit procedure carried out in the frequency domain. Such 
a process can produce a set of discrete autoregressive moving average (ARMA) filters which 
together reproduce the multichannel system frequency response to the given accuracy. The 
fitting process is based on a modified version of the vector fitting algorithm [T7] adapted 
to work in Z-domain [18] . This procedure allows the frequency response of the coloring filter 
to be fit with ARMA functions expanded in partial fractions: 



H(a;) 



hu (w) 



hiM (w) 



^H(z) 



hMl (w) ■ ■ ■ hMM (w) y 

N 

^'ij,k 



( hu (z) ■■■ hiM (z) ^ 
y hMl (z) ■■■ hMM (z) J 



k=l 



Pij,kZ 



(18) 



The number of poles required to obtain a satisfactory fit of the model transfer function is 
automatically determined by an iterative procedure in which the number of poles is increased 
by one at each step of the fit loop. The iteration stops when the mean square error between 
fit function response and model response comes to a value smaller than the user-defined 
threshold. Since the fit is performed in the Z-domain, the noise generation procedure turn 
out to be free from aliasing as discussed in section |T} 



It is worth noting that the eigenvectors in equation (16) are defined up to an arbitrary 
phase factor. This means that the columns of the V {u) matrix can be multiplied by an 
arbitrary phase factors e*"^ without changing their property of being eigenvectors of the 
cross-spectral matrix. Such phase arbitrariness does not extend to the single elements in the 
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columns of V (u); they are elements of the same eigenvector of V (u) and their phase relation 
must be carefully preserved during the fit process because it is connected to the correlation 
properties of the multichannel system. It can happen, after the eigen-decomposition process, 
that the phase of the elements of V (u) is such that the frequency response of the elements of 
H (cj) cannot be fit with stable poles. In that case, the poles must be stabilized after the fit 
process by the application of an all-pass filter [TT] . The all-pass function substitutes unstable 
poles with the inverse of their conjugates which are stable. Its magnitude (absolute value) is 
1 at each frequency. As already mentioned, the phase relation between the elements of each 
column of H [z) must be kept constant to prevent the corruption of the system correlation 
properties. The proper all-pass filter for the elements of H {z) stabilizes the unstable poles 
of the given hij (z) and at the same time adds an extra phase in order to keep the phase 
relation between the columns of H {z) constant. Therefore each element of H {z) is modified 
according to: 




(19) 



Here, NJ^a is the number of unstable poles in has (z). The function Ff ( ^ ) 
has the purpose of poles stabilization for the element hap (z) of H(2;). The product 



n 

for t 



-yfi 



5i v^p;^,.-i 



provides an extra phase coming from the poles stabilization procedure 

le other elements of the same column of H {z). In this way, each element of the matrix 
H {z) comes with such a phase allowing the representation with just stable poles and, at 
the same time, the original phase relation between elements of the same column of H {z) is 
preserved. A second fit step with stable poles provides usable filters for noise generation. 



C. Filter Initialization 



The output of a linear discrete causal filter is a function of the present and all previous 
input values. Since the input to the filter must start at some time t = 0, the output process 
will consist of an unwanted filter transient response added to the desired stationary random 
process. One possible approach to handle the problem of the filter transients is to wait for 
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FIG. 1: Scheme for the apphcation of a multichannel filter expanded in partial fractions. 



the time necessary for the transients to decay to an acceptable level. However the transient 
response is proportional to the filter impulse response and if there are poles too near the 
unitary circle of the complex plane the transient response could last for an unacceptably 
long time. Moreover the practice of hand removing initial data is always inaccurate and 
can hide important features or introduce fake signals especially when the spectrum spans 
several decades in frequency. Therefore the filter for the noise generation should be always 
properly initialized. 

We are searching for an initialization for the recursive equation of a ARMA process 
written as sum of partial fractions (Figure ([T]), Equation (18)): 



Xij,k (n) = Pij,kXij,k {n-1) + Tij^kSj [n) 
Vijin) = ^Xij^kin). (20) 



k=l 



The recursive equation (20) calculates the output of the system at a given time on the 
basis of the input ej (n) at the same time and the information from the previous output 
Xij,k {n — 1). In the present case the input process ej {n) is an element of a discrete multi- 
channel unitary variance white noise process such that: 
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(21) 



In order to properly initialize the recursive equations for the implementation of the mul- 
tichannel filter, it is necessary to calculate the covariance matrix of the initial states Xij^k (0). 



It should be considered that the M processes represented in equation (14) are not inde 



pendent from each other. A combined process should be defined in which all the recursive 
equations (equation 20) for each of the processes are incorporated. Readily it is seen 



that, thanks to the delta correlation properties (21) of the input signals, the processes ap 



plied to different input data series are independent. This means that, instead of building 
a single combined process, one has to build M independent processes which combine those 
applied to the same input. The new processes can then be written as: 



Xij,Nij (n) = pij^Ni.Xij^Nij {n-1) + rij^N^^Ej (n) 



XMj,i {n) = PMj,iXMj,i {n-1) + rMj,iej {n) 



, XMj,NMj (n) = PMj,NM,XMj,N,,, {n-l)+ rMj,NMj£j i^) 



(22) 



The covariance of such processes (22) can be written as: 



{Xi,a (n) {m)) = {pi,aXi,a {n - l)p*„Xj,/3 ("^ - 1)> + {ri,ae^ (n) rl^e* (m)) 

a = l,...,Nu + --- + NMi 

P = 1, . . . , A^i, + ■ ■ ■ + Nmj 

2,j = l,---,M. (23) 

Where again, the symbol () represents the expectation value operator. Assuming sta- 
tionary processes: 

{Xi,a (n) xl)3 ("^)> = Rij,af3 (n, m) = Rij^^ii {n - m) , (24) 
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and thanks to the properties of the input functions (21): 



R 



ij,al3 



[n — m] 



3,1 



-6ij6 {n — m) , 



1 - Pi,aPj,/3 

which provides the desired covariance for the first state of the recurrence sequences: 



(25) 



Rj,af5 (0) 



(26) 



Initial states for the filter recurrence sequence can then be generated by a multivariate 



noise generator according to the given covariance (26). As an alternative, they can be 



calculated from random independent variables through a linear transformation [6j of the 
type: 



Xj (0) = A,- ■ 77, 



J' 



(27) 



where r]j is a column vector of Nij + ■ ■ ■ + N^j independent zero mean and unit variance 
random numbers and Aj is a {Nij + ■ ■ ■ + Nmj) x {Nij + ■ ■ ■ + Nmj) transformation matrix. 



If Rj,ai3 (0) is calculated for the variables in equation (27), it is found: 



R, (0) = (A, ■ 77,) ■ (A, ■ ri.y 

= A, ■ I ■ At, (28) 

and it is readily seen that: 

A, = V, • (29) 

where and Sj are the eigenvector and eigenvalue matrices of (0). 

III. A CASE STUDY, LTP ALONG X AXIS 
A. Response model and fit 

An application of the noise generation procedure is presented for a two channel system 
simulating the LISA Technology Package (LTP) along the principal measurement axis [T]- 
[T3HT5] . The complete set of algorithms are available as MATLAB tools in the framework 
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of the LTPDA toolbox [TU [19] and can be freely downloaded, together with the complete 
toolbox, at the LTPDA project web page pU] . 

The expected power spectra and cross-power spectrum at the output of the system can 
be calculated (Figure (|2])) on the basis of some assumptions on the properties of input noise 
sources |4j. The noise coloring filters can be calculated following the procedure described 



in paragraphs II A and II B A frequency domain fit is performed on the models for the 
coloring filters obtained by eigendecomposition (frequency by frequency) of the expected 
cross-spectral density matrix. The fit procedure takes around 200 seconds on a standard 
desktop machine [21]. The four transfer functions hn {z), hi2 (z), /i2i {z) and {z) have 
respectively 25, 30, 28 and 30 poles. The fit loop stops when the mean square error between 
fit function response and model response is smaller than 1 x 10~^. The response of the filter 



designed to reproduce the cross-spectral density can be calculated according to equation ( 15 ). 
It can then be compared with the model cross-spectral density of the system as reported in 
figure [3] 
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FIG. 2: Model power spectra and cross-spectrum for the signals at the output of the two 
channels. 5*11 and 5*22 are real values so they do not appear in the bottom plot. 



In order to compare the correlation properties of expected model with the fit results, it 
is useful to introduce the complex cross-coherence: 
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FIG. 3: Comparison between model power spectral densities and fit result, a) Output of 
the first channel, b) Output of the second channel. 



5*12 (uj) 



(30) 



^/Sii (w) S22 {^) 

where 5*12 (w), 5*11 (w) and S22 (w) are cross-spectrum and power spectra of the first and 
second channels. The real and imaginary part of the cross-coherence for the expected and 
fit cross-spectral matrices are reported in figure |4} 
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FIG. 4: Comparison between expected coherence and fit result, a) Real part, b) Imaginary 

part. 



Any discrepancy between the expected model and the fit model can be considered as 
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a systematic error in the procedure, whose influence on the process can be minimized by 
increasing the fit accuracy. Clearly this has a computational cost in terms of the number of 
poles required to match the accuracy goal and on the amount of time required to complete 
the fit loop. Hereafter, the fit model will be considered as the reference model. 



B. Noise generation tests 

Once the two channel noise coloring filter is obtained it can be used to generate a two 



channel noise data series according to the procedure described in paragraph |II C[ Data series 
are 3 x 10^ seconds long at a sampling rate of 10 Hz. The chosen rate is the same as the 
LTP experiment operations, e.g., the control forces acting on test masses and the spacecraft 
will be calculated by controllers on the basis of 10 Hz sampled data streams. 

In order to realize a statistically meaningful test, = 500 independent realizations of the 
two channel process were generated. Power spectra of the two channels are calculated with 
the windowed periodogram method using a 4-term Blackman-Harris window [21] . The choice 
of such a window is justified by the requirements in terms of spectral leakage performances. 
A 4-term Blackman-Harris window, having the highest side-lobe level of -92 dB (relative to 
the main lobe level) [21], is indeed one of the best-performing available window in terms of 
spectral leakage suppression. The realizations of the power spectrum were averaged and 
compared with the reference model expectation (figure [s]). 



The reported uncertainty is calculated under the assumption that (Tmean (w) 



-'pop 



{UJ) 



where apop (cu) is the sample standard deviation of the spectra population at a given fre- 
quency. In doing this we have considered that, since the power spectrum of a X2 distributed 
variable |22j, the mean of A^ independent realization of the same spectrum will also be X2N 
distributed. Since in our case A^ = 500, the average of the spectra is distributed with 1000 
degrees of freedom, and such a variable can be considered to be Gaussian distributed with 
reasonable accuracy [25]. The expected standard deviation (normalized to the mean) for the 
equivalent normal distribution is 0.045 where we measure on average amean = 0.044 ± 0.002 



In the procedure for the calculation of the power spectrum, data are multiplied in the 
time domain for the time response of the window function. As this operation corresponds 
to a convolution in the frequency domain, the reference model must include also the effect 
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FIG. 5: Averaged power spectral density compared with reference model, a) First channel. 

b) Second channel. 



of the window function. Windowed spectra can be calculated as: 



27Tk 



2nN 



s{n) 



9=0 



(31) 



where is the number of samples in the data series and Wq are the time samples of the 



window function. The integral in equation (31) is numerically evaluated, and the results are 



reported in figure [5j As can be seen, the effect of the window convolution is visible at the 
lowest frequencies, were the departure from the reference model is remarkable. 

A quantitative analysis of the results is better performed with the introduction of the 
variable: 



AS(u) 



Syy {Uj) - S {U) 



(32) 
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S (w) represents the expected value for the spectrum (at each frequency) and Syy (u) is 
the estimated spectrum (averaged over 500 reahzations) at the given frequency. AS can be 
considered distributed in accordance to a function. Therefore its expectation value 

is 0. 

AS is calculated for the simulated data and the windowed model with respect to the 
reference model. Results are reported in figure |6j A 99.97% confidence interval is calculated 
on the basis of the statistical properties of AS*. 

The effect of the window on the spectra calculation is noticeably high on the first 3 
frequency bins, where the deviation from the reference model exceeds the confidence interval. 
In addition, it is clearly observable up to the 10*'^ bin. The windowed model and the 
simulated data are consistent on the basis of the chosen confidence region. 

A considerable number of data points, especially at high frequencies, lie outside the 
confidence levels, and it is of fundamental importance to assess if such outliers are caused 
by the random nature of the data, or if they come from systematic errors in the spectral 
estimation or data generation processes. The confidence levels at 99.97% define a region in 
which the data are expected to lie with that probability. This also means that in 0.03% of 
the observations an outlier can be observed. As we are dealing with datasets of 1.5 x 10^ 
points the number of expected outliers is high. 

In order to distinguish between systematic outliers and statistical outliers the averaging 
process over 500 independent realizations was repeated 5 times and the frequencies at which 
the values of AS were outside the defined confidence interval were recorded. The first 10 
frequency bins are excluded from the numbering because they are systematically affected by 
the spectral widow effect. Figure [7] reports a histogram of the cumulative count of the outliers 
frequencies for the 5 different realizations. If an outlier is originated by a systematic error, 
then it is expected to be counted 5 times. As can be observed from figure [7], the maximum 
value obtained is 2 for both channels; this is a definitive indication of the statistical nature 
of the observed outliers. 

As stated above, the correlation properties of the data series can be explored with the 



sample coherence calculated as in equation (30). Power spectra and cross-spectrum were 
estimated with the averaged Welch periodogram method using a 4-term Blackman-Harris 
window over 145 data segments 6 x 10^ points long. The separate 500 realizations are then 
averaged and, assuming the averaged process is approximately normally distributed, the 
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FIG. 6: AS* [u) calculated for a) First channel and b) Second channel. Simulated results 
are compared with the model expectation and the windowed model expectation. A 99.97% 
confidence interval is calculated for quantitative comparison with the models. 

error on the estimation is calculated as described above. Results are reported in figure 
[8] and compared with the expectation from the reference model. Since the coherence is 
constructed from a ratio between cross-spectrum and power spectra the effect of the window 
on the lowermost frequency bins is strongly attenuated. 

The simulated data (averaged over 500 realizations) and the reference model are in sat- 
isfactory agreement within the tolerance region defined by the uncertainty. On the basis 
of the above discussion, the oscillations observed in the coherence curves (figure [s]) can be 
associated with the statistical fiuctuations caused by the random nature of the data. 
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FIG. 7: Histogram of the outliers frequencies for ASi and AS2. 
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(a) (b) 

FIG. 8: Averaged cross-coherence compared with the reference model, a) Real part, b) 

Imaginary part. 

IV. CONCLUSIONS 

A robust procedure for the generation of multichannel stationary noise with a given cross- 
spectral matrix is reported. Based on some assumptions on the noise sources acting on the 
system under study, an expected model for the cross-spectral matrix of the multichannel 
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output noise can be developed. Prom such a model the noise coloring filters are identified by 
an eigendecomposition of the cross-spectral matrix (frequency by frequency) and a frequency 
domain fit procedure. A multichannel colored noise data series can then be generated from 
a multichannel 6 correlated random noise process provided that the recurrence equations are 
properly initialized in order to avoid transients at the beginning of the noise sequence. It is 
demonstrated that the only source of systematic errors in the process is associated with the fit 
procedure. On the other hand, the accuracy of the fit can be increased at the expense of the 
computational cost of the whole process; this, in principle, ensures that the process reaches 
the desired accuracy. An average over 500 independent realizations of the multichannel 
noise process has demonstrated the statistical consistency between generated noise and the 
reference model if the effect of the spectral window is taken into account. Oscillations 
in the averaged spectra with respect to the model can be unambiguously attributed to 
statistical fiuctuations. The analysis reported demonstrates that the tool can be applied 
for the calibration of spectral estimators in experiments where noise spectral energy content 
must be estimated with very high accuracy, as is the case for the LTP experiment. MATLAB 
based algorithms are available for free download at the LTPDA project web page. 
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